library(tidyverse)
theme_set(theme_classic())
knitr::opts_chunk$set(warning = F, message = F)

cargo datos

Transformo la variable id en un factor porque no es numérica realmente. Transformo las variables character en factor manualmente, ya que readxl::read_excel no tiene argumento stringAsFactors.

Estas transformaciones las hago porque tiene más sentido la descriptiva sobre factores. En todos los casos las variables pueden tomar uno de un conjunto finito de niveles.

d <- readxl::read_excel("dataset.xls",) %>% janitor::clean_names() %>%
  mutate(across(where(is.character), factor),
         id = factor(id),
         color = factor(color, levels = c("ROJO", "VERDE", "AZUL", "VIOLETA")))

Tipos de variables

Categóricas: id, color, tiene_pin, vt, cliente_pas, churn

Cuantitativas: edad, antiguedad, consumo_tc, mov90_cta, sueldo

descriptiva

En las variables categóricas vemos que no tenemos ningún NA, lo cual es bueno. Sólo tenemos un problema en la variable tiene_pin, que debería ser SI/NO, pero hay una tercera opción, “2”, con 65 observaciones. El resto de las variables tienen buena pinta a orden cero.

En la variables cuantitativas, la variable edad tiene al menos un outlier, de 900 años, y lo mismo para la variable antiguedad, con 337. Luego en la variable consumo_tc tenemos al menos un outlier: un consumo del orden de 10^9, lo cual probablemente no tenga sentido. La columna mov90_cta tiene un máximo bastante mas alto que la media, pero en principio podría ser una cuenta muy utilizada (comercial, por ejemplo). Por último, en la columna sueldos tenemos muchas observaciones faltantes, con 13298 NA.

d %>% skimr::skim()
Data summary
Name Piped data
Number of rows 18179
Number of columns 11
_______________________
Column type frequency:
factor 6
numeric 5
________________________
Group variables None

Variable type: factor

skim_variable n_missing complete_rate ordered n_unique top_counts
id 0 1 FALSE 18179 368: 1, 398: 1, 406: 1, 511: 1
color 0 1 FALSE 4 VER: 7371, AZU: 6361, ROJ: 2947, VIO: 1500
tiene_pin 0 1 FALSE 3 SI: 13498, NO: 4616, 2: 65
vt 0 1 FALSE 2 NO: 11365, SI: 6814
cliente_pas 0 1 FALSE 2 NO: 13298, SI: 4881
churn 0 1 FALSE 2 NO: 17987, SI: 192

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
edad 0 1.00 43.56 17.41 12 31 41 54.00 900 ▇▁▁▁▁
antiguedad 0 1.00 69.95 73.66 1 10 39 111.00 337 ▇▃▂▁▁
consumo_tc 0 1.00 21134456.71 327042640.71 0 0 0 4471.59 9509228475 ▇▁▁▁▁
mov90_cta 0 1.00 43.72 57.14 0 3 20 66.00 826 ▇▁▁▁▁
sueldo 13298 0.27 31877.97 58745.72 10 14047 22553 37471.00 2480722 ▇▁▁▁▁

biplots

Primero hago un gráfico con todas las combinaciones, y coloreando según el la variable color. Esta figura, asi como está, no es muy informativa, ya que es muy difícil de interpretar, por lo que más abajo la vamos a ir partiendo.

d %>%
  select(-id) %>%
  GGally::ggpairs(aes(color = color), showStrips = T)

cuantitativas

Por un lado graficamos las variables cuantitativas:

d %>%
  select(where(is.numeric)) %>%
  GGally::ggpairs()

En la columna correcpondiente a edad vemos claramente un grupo de observaciones muy extremas. Podemos graficarlo como boxplot, separando por color, para facilitar la visualización. Vemos que hay outliers de edades en todas las categorías.

Asumiendo que se tratan de personas físicas, y dado que la persona mas longeva en Argentina tiene apenas 114 años, usamos ese número como cutoff, y eliminamos las observaciones superiores. Probablemente sea una sobreestimación, pero prefiero una aproximación conservadora.

d %>% ggplot(aes(color, edad, color = color)) +
  geom_boxplot() +
  scale_y_continuous(n.breaks = 10)

d %>% filter(edad <= 114) -> d

Grafico las cuantitativas nuevamente y vemos que la columna de edad está mas linda.

d %>%
  select(where(is.numeric)) %>%
  GGally::ggpairs()

Con respecto a antiguedad, vemos que hay un subgrupo muy separado de los otros (por sobre los 300 meses). Grafico en boxplot, separando por color y graficando los outliers independientemente, con un poco de jitter, para tener una idea de que pasa.

En las cuatro clases hay una concentracion de antiguedades alrededor de los 284 meses y otra a los 337, lo cual de mínima es llamativo. Un hipótesis es que tengamos un faltante de datos con antiguedades entre 284 y 337. Hipotetizando más aún, quizá esas antiguedades intermedias fueron colapsadas en 284 y 337.

d %>% ggplot(aes(color, antiguedad, color = color)) +
  geom_point(aes(alpha = I(antiguedad) / 700), position = "jitter", color = "gray") +
  geom_boxplot(outlier.shape = NA, alpha = 0) +
  geom_hline(yintercept = 284, linetype = "dashed", alpha = 0.3) +
  geom_hline(yintercept = 337, linetype = "dashed", alpha = 0.3) +
  scale_y_continuous(n.breaks = 10)

Con respecto a la variable consumo_tc, aproximadamente la mitad de los clientes en la base no usaron la tarjeta en el mes de análisis (o al menos no está reportado). Para el resto, transformo logaritmicamente para facilitar la visualizacion de las distribuciones en un boxplot, de nuevo separado por colores. Ademas, grafico los datos encima. Vemos que el alto consumo de tarjeta no era un outlier! Quizá los numeros enormes puedan explicarse por la unidad usada, que desconocemos.

Otra opción es que todas esas observaciones sean outliers, por error multiplicadas por un factor grande.

filter(d, consumo_tc > 0) %>% ggplot(aes(color, consumo_tc)) +
  geom_point(aes(alpha = I(log(consumo_tc)) / 50), position = "jitter", color = "gray") +
  geom_boxplot(alpha = 0) +
  scale_y_log10()

Con respecto a las variables sueldo y mov90_cta si bien hay outliers definidos “a la Tukey”, ninguno de los valores son “imposibles” por lo cual no hago nada al respecto.

d %>% ggplot(aes(color, sueldo)) +
  geom_boxplot()

d %>% ggplot(aes(color, mov90_cta)) +
  geom_boxplot()

categoricas

Como mencionamos antes, la variable tiene_pin tiene un nivel extra incorrecto. Lo elimino para simplificar el gráfico (son sólo 65 observaciones, en 18000).

d %>%
  select(where(is.factor), -id) %>% 
  filter(tiene_pin != 2) %>% 
  GGally::ggpairs(aes(color = color), showStrips = T)

Graficando ese grupo por separado, contra todas las cuantitativas, vemos que el unico detalle notable es que todos los casos corresponden a cuentas cuyos sueldos no están en la base de datos. Pero notable hasta ahí, porque esa columna tiene un 26% de complete_rate.

d %>%
  select(where(is.numeric), tiene_pin) %>%
  filter(tiene_pin == 2) %>%
  GGally::ggpairs()

El resto de las categóricas no parece tener nada raro. Y al graficarlas contra las cuantitativas nada me llama la atención particularmente.

outliers

En resumen, tenemos outliers en la columna edad, que eliminamos tomando como cutoff la edad de la persona mas longeva alguna vez registrada en Argentina.

En la columna tiene_pin tenemos unas ~65 observaciones con un nivel mal asignado. El resto de las variables en esas observaciones tienen pinta normal, por lo que les asignaria NA en la variable tiene_pin y listo.

Por último, me resultan extrañas las observaciones de la variable consumo_tc; son montos demasiado altos. Ademas, no correlacionan con la variable sueldo lo cual es extraño.

d %>%
  filter(consumo_tc < 1e6) %>%
  select(sueldo, consumo_tc) %>%
  GGally::ggpairs()

Graficando un scatterplot, con los ejes logaritmicos vemos que el grupo de consumo_tc alto es tiene forma similar al del resto, pero como “subido”. Lo cual soporta la hipótesis de que son todos outliers, y que por error les agregaron algo, o los multiplicaron por algún factor. (nótese que no hay puntos rojos en la nube inferior, ya que el orden de las capas de ggplot hace que queden por encima de los puntos negros.)

d %>% ggplot(aes(sueldo, consumo_tc)) +
  geom_point() +
  geom_point(data = filter(d, consumo_tc > 1e6), color = "red")

d %>% ggplot(aes(sueldo, consumo_tc)) +
  geom_point() +
  geom_point(data = filter(d, consumo_tc > 1e6), color = "red") +
  scale_y_log10() +
  scale_x_log10()

Con respecto a estos datos, quizá les asignaria NA en la variable consumo_tc.

mahalanobis

Para calcular las distancias de Mahalanobis necesito unicamente las variables cuantitativas. Ademas, tengo que eliminar los NAs.

d %>%
  select(where(is.numeric)) %>%
  filter(complete.cases(d)) %>% 
  mahalanobis(center = colMeans(.), cov = cov(.)) -> mahalanobis_distances

d %>% 
  filter(complete.cases(d)) %>%
  mutate(mahalanobis = mahalanobis_distances) -> d_distances

Graficamos unos lollipops con las distancias, y tomando el indice de la observacion como x. Vemos que hay observaciones con distancias muy grandes.

d_distances %>%
  mutate(idu = as.numeric(rownames(d_distances))) %>% 
  ggplot(aes(idu, mahalanobis)) +
  geom_point() +
  geom_segment(aes(xend = idu, yend = mahalanobis, y = 0))

Para definir un cutoff, puedo usar la distribución de distancias, y tomar como outlier a las observaciones con distancias en el 5% mas alto de la distribución.

d_distances %>%
  ggplot(aes(mahalanobis)) +
  geom_density() +
  geom_vline(aes(xintercept = quantile(mahalanobis, 0.95)), linetype = "dashed") +
  scale_x_log10()

Graficamos nuevamente los lollipops, pero agregando el cutoff. Grafico el eje y logaritmico, ya que el rango es muy grande, y no se aprecia la posicion del cutoff. Elimino las lineas del lollipop para evitar el overplotting.

d_distances %>%
  mutate(idu = as.numeric(rownames(d_distances))) %>%
  ggplot(aes(idu, mahalanobis)) +
  geom_point(alpha = 0.5, size = 1) +
  scale_y_log10() +
  geom_hline(aes(yintercept = quantile(mahalanobis, 0.95)),
             linetype = "dashed",
             size = 1, color = "red")

Finalmente, grafico el consumo de tarjeta de credito en funcion del sueldo, y coloreo los outliers en rojo. Vemos que la mayoria de los puntos de la nube superior son efectivamente outliers. En la nube inferior tambien hay outliers, probablemente debido a las distancias con respecto a las otras variables.

d_distances %>% ggplot(aes(sueldo, consumo_tc)) +
  geom_point() +
  geom_point(data = filter(d_distances, mahalanobis > 10), color = "red") +
  scale_y_log10() +
  scale_x_log10()

correlogram

Para las variables cuantitativas podemos hacer un correlograma

d %>% 
  select(where(is.numeric)) %>% 
  GGally::ggcorr(label = T, hjust = 0.7, vjust = 2, angle = -45) +
  theme_minimal()

profile plots

Para los graficos de perfiles, me interesan las medias/medianas para diferentes agrupamientos segun las categoricas.

colors <- c("red", "green3", "blue", "violet", "black", "gray")

d %>%
  select(-id) %>%
  filter(tiene_pin != 2) %>%
  pivot_longer(
    cols = c(color, tiene_pin, cliente_pas, churn),
    names_to = "g_var", values_to = "g_val"
  ) %>%
  
  group_by(g_var, g_val) %>%
  summarise(across(is.numeric, mean, na.rm = T)) %>%
  pivot_longer(-c(g_var, g_val)) %>%
  
  ggplot(aes(name, value, color = g_val, group = g_val)) +
  geom_point() +
  geom_line() +
  scale_color_manual(values = colors) +
  facet_wrap(~g_var, scales = "free_y") +
  labs(x = "") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

colors <- c("red", "green3", "blue", "violet", "black", "gray")

d %>%
  select(-id) %>%
  filter(tiene_pin != 2) %>%
  pivot_longer(
    cols = c(color, tiene_pin, cliente_pas, churn),
    names_to = "g_var", values_to = "g_val"
  ) %>%
  
  group_by(g_var, g_val) %>%
  summarise(across(is.numeric, median, na.rm = T)) %>%
  pivot_longer(-c(g_var, g_val)) %>%
  
  ggplot(aes(name, value, color = g_val, group = g_val)) +
  geom_point() +
  geom_line() +
  scale_color_manual(values = colors) +
  facet_wrap(~g_var, scales = "free")

sample

Primero calculo la proporcion de churn en todo el dataset.

set.seed(1)
d %>%
  select(churn) %>%
  table() %>%
  prop.table()
## .
##         NO         SI 
## 0.98943488 0.01056512

Luego tomo una muestra del 20%. Vemos que da practicamente igual.

d %>% 
  select(churn) %>% 
  sample_frac(0.2) %>% 
  table() %>%
  prop.table()
## .
##          NO          SI 
## 0.990096286 0.009903714